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ABSTRACT 

We present a maximum likelihood method for determining the spatial proper- 
ties of tidal debris and of the Galactic spheroid. With this method we characterize 
Sagittarius debris using stars with the colors of blue F turnoff stars in SDSS stripe 
82. The debris is located at {a, 5,R) = (31.37° ± 0.26°, 0.0°, 29.22 ± 0.20 kpc), 
with a (spatial) direction given by the unit vector < —0.991 ± 0.007 kpc, 0.042 ± 
0.033 kpc, 0.127±0.046 kpc >, in Galactocentric Cartesian coordinates, and with 
FWHM = 6.74 ±0.06 kpc. This 2.5°-wide stripe contains 0.9% as many F turnoff 
stars as the current Sagittarius dwarf galaxy. Over small spatial extent, the debris 
is modeled as a cylinder with a density that falls off as a Gaussian with distance 
from the axis, while the smooth component of the spheroid is modeled with a 
Hernquist profile. We assume that the absolute magnitude of F turnoff stars is 
distributed as a Gaussian, which is an improvement over previous methods which 
fixed the absolute magnitude at M go = 4.2. The effectiveness and correctness of 
the algorithm is demonstrated on a simulated set of F turnoff stars created to 
mimic SDSS stripe 82 data, which shows that we have a much greater accuracy 
than previous studies. Our algorithm can be applied to divide the stellar data 
into two catalogs: one which fits the stream density profile and one with the 
characteristics of the spheroid. This allows us to effectively separate tidal de- 
bris from the spheroid population, both facilitating the study of the tidal stream 
dynamics and providing a test of whether a smooth spheroidal population exists. 
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Introduction 



1.1. The Milky Way's Spheroid 



In the late 1980s, the Milky Way's spheroid population was described as slowly rotat- 
ing, with density distribution given by p oc r -3,5 ( IFreemanl 1 198 71 ) . There was a controversy 
concerning the flattenin g of the spheroid, since kinematic s tudies of the spheroid stars sug- 
gested a flatter spheroid (IGilmore. Wvse. and Kuijkenlll989l ) and star counts often suggested 
a more spherical spheroid (Bahcall 19861 ). By early in the 21 st century, the spheroid was still 
thought of as a smooth power law distribution, but studies were starting to show that the 
shape of the spheroid depended on the type of star being observed, and it was noted that 
at least some of the spheroid, if not all, was comp osed of debris from hierarchical struc- 
ture formation (iFreeman fe Bland-Hawthornd 120021 ). The Sagittari us (Sgr) dwarf galaxy 
( Ibata. Gilmore. and Irwrr] ~ 1994f) and its associ ated tidal stream ( Yanny. Newberg et al. 



200( 



o|; Ibata. Lewis, et al 



2001 



Ibata et al.l 120011 ) was the one known example of merging 
in the present day. In the last ten years, the discovery of substructure has dominated the 
discussion of the Galactic spheroid. The discovery of substructure has been driven primarily 
by the Sloan Digital Sky Survey (SDSS) and the related Sloan Extension for Galactic Under- 
standing and Exploration (SEGUE), but other surveys such as the Quasar Equatorial Survey 
Team (QUEST) and the Two Micron All Sky Survey (2MASS) have also been influential. 
Discoveries include new tidal debris streams, dwarf galaxies, and globular clusters. 

The positions, velocities, and met allicities of Sgr dwarf spheroidal tidal debris stars 



have been measured all across the sky (INewberg. Yanny et al 



20021; 



Bellazzini et al. 



20031: 



20041 ; iMartinez-Delgado et al. 



Maiewski et al.ll2003l; INewberg. Yanny et al.ll2003l ; iMajewski et al.l . 

2004l : lFellhauer et aliEoodlChou et al.ll2007h. These measurements have th e n been compared 



with models of tidal disruption Jjohnston et al" 



1998 



2004 



Gomez-Flechoso et al. 1999; Helmi and White 



1995 



Ibata et al. 



2001 



1997 



Johnston et al 



Ibata and Lewis 



2002: 



Law et al 



2005[ l. Given that the data spans three spatial dimensions plus radial velocities, but 



both the models and data are presented in papers primarily as two dimensional plots, it 
has been difficult to make detailed comparisons in multidimensional space. In this paper 
we develop a technique that would allow us to create a catalog of stars that has the same 
spatial distribution as the stars in the Sgr tidal stream. That catalog could then be made 
available to the modeling community, who would be able to transform the coordinates into 
whatever system they find most natural to make comparisons. 

At least three newly discovered tidal debris streams are thought to be associated with 



dwarf galaxies, including: the Mon oceros stream in the Galact ic plane (INewberg. Yanny et al. 



20021 : lYannv. Newberg et al.ll2003 ). the Virgo Stellar Stream flVivas et al.ll200ll : iDuffau et al. 
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2006 



Newberg et al. I l2007h . and the Orphan Stream (Grillmair 2006) in the Field of Streams (IBelokuroy et 



2006bJ). An additional piece of tidal debris was found in Triangulum- Andromeda by lRocha-Pinto et al 



( 120041 ) which might be part of the Monoceros stream (IPennarrubia et al.l 120051). Tidal tails 



spanning many tens of degrees across the sky have bee n found around Pa l 5 (IQdenkirchen et al 



2001 



2006 



Rockosi et al. 



2002 



Gril 



Belokurov et all2006d ). 



mair & Dionatos 2006), and NGC 5466 ( Grillmair fc Johnson 



Grillmair fc Dionatosl (120061 ) find a similar-looking tidal stream 

jright 



for which there is no known progenitor globular cluster. Eight ne w low surface 



ness dwarf galaxy satellites of the Milky Way have been discovered (IWillman et al. 



Zuckcr et al. 2006a: Be 



2007a 



Irwin et al 



okurov et al1l2006al : IZucker et al1l2006bl ; lGrillmairil2006l 



2005a 



Belokurov et al 



20071 ). nearly doubling the number of known Milky Way dwarfs. IBelokuroy et al 



( I2007bl ) suggest that there is a new structural component of the spheroid called the Hercules- 



Aquila cloud. 

With all of this substructure detected in the densities of spheroid stars, one wonders 
how much of the spheroid is smooth, and what the shape of the smooth spheroid might be. It 



is no wonder that standard Galactic models, such as the the Besangon model (jRobin et al. 



20031 ). which uses a power law density pro file, do not m atch the SDSS star counts. The 
spheroid s hape is a current con troversy as iHelmil (120041 ) has found prolate models favor- 
able while Johnston et al. (l2005h has found oblate models best. The large number of very 
significant over-densities of stars in the spheroid show that pencil beam studies and stud- 
ies that extrapolate from the solar neighborhood could be contaminated by local structure 
and may not be indicative of the global shape of the Milky Way spheroid. Recent results 
from F turnoff stars from the SDSS suggest the smooth component of the spheroid might 
be asy mmetric about the Galactic center, so that no axial l y symmetric spheroid model can 
befit JNewberg and Yannvi [20051 . bood : Isavage et alJliooS ku.Deng fc Hull2006h . However, 
accurate measurements of the smooth component require that spatial substructure be iden- 
tified and quantified — otherwise it is difficult to know which stars to fit to smooth spheroid 
models. 

We need to develop a model that is as complex as the data. In this paper we start this 
process by developing a maximum likelihood method that can be used to fit a smooth stellar 
spheroid containing one debris stream, using one 2.5° stripe of SDSS data, and show that our 
algorithm produces the correct results on a simulated dataset. In the future, we will extend 
the algorithm to include larger datasets with multiple pieces of debris. The long-term goal 
is to develop a model of the spheroid that fits the observed, lumpy density distribution of 
stars in the spheroid. 
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1.2. The Sloan Digital Sky Survey 



The SDSS is a large, international collaboration that was originally established to find 
the largest structures of galaxies in the Universe from an imaging survey of 10,000 sq. deg. of 
sky and spectra of 1,000,000 galaxies selected from the photometry (from images). The SDSS 
can simultaneously obtain 640 spectra per pointing, using two 320-fiber double (blue/red) 
spectrographs. The angular position of a galaxy in the sky is easy to measure from the 
imaging survey, and for distant galaxies the distance is well estimated from a spectroscopic 
measurement of the Hubble redshift. By combining the imaging with the spectroscopic 
survey, the SDSS has built up a three dimensional picture of the distribution of galaxies 
in the U niverse using a d edicated 2.5 meter telescope at Apache Point Observatory in New 
Mexico jYork et al.ll200oh . 



Images are obtained with an array of thirty 2048 x 2048 pixel CCD cameras operated 
in a drift-scan (time delay and integrate, TDI) mode, which produces six "scan lines" of 
imaging data that are each 13.6' wide and grow longer at a rate of 15 degrees per hour for 
the length of that "run." In each of the six "scan lines," the sky is imaged in five optical 
filters: u, g, r, i, and z; the time at which each astronomical object is imaged is a few minutes 
different for each passband. When a second set of six scan lines is observed that fill in the 
gaps between the first set of six scan lines, they produce a contiguous "stripe" of data of 
width 2.5° and length that depends on the length of time the sky was observed in the "runs." 
Each 2.5°-wide stripe follows a defined great circle on the sky, and is built up of two or more 
"runs" of data, as many runs as are needed to traverse the SDSS survey area. 

The entire sky is divided up into 144 numbered stripes that start and end at the survey 
poles: (l,b) = (209.33°, -7°) and (I, b) = (29.33°, 7°). Stripes 10 and 82 are centered on 
the Celestial Equator, with stripe 10 in the North Galactic Cap and stripe 82 in the South 
Galactic Cap. The other stripes are sequentially numbered with inclinations 2.5 degrees 
apart. If the entire length (180°) of every stripe were imaged, we would have 64, 800 sq. 
degrees of imaging data for a sky that is 41,253 sq. deg. The overlaps between stripes 
increase towards the survey poles. Since the SDSS generally observes only parts of the sky 
that are more than 30° from the survey poles, there is about 50% overlap at the ends of 
the stripes and very li ttle overlap on the survey equat or, which is at RA= 185°. With the 
release of SDSS DR6 (lAdelman-McCarthy et al.l 120071 ). a contiguous 8500 sq. deg. area of 
the North Galactic Cap is now publicly available. It is pieced together with 31 adjacent SDSS 
stripes. Three stripes of data are available in the South Galactic Cap, including stripe 82 
on the Celestial E quator. Further i nform at ion about the SPSS, in cl uding survey ge o metry , 
can be found in: IStoughton et al.l (|2001l). lAbazaiian et al.l (120031). iFukugita et al.l (119961 ). 



Gunn et al.l (Il998f ). iHogg et al.l (I200lh . IPier et all (12003 ) 



Smith et al 



(j2002h . 
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Because the imaging survey is large and well calibrated, and because the spectroscopic 
survey included many Galactic stars, in some cases by design and in some cases accidentally, 
the SDSS has made a significant contribution to our knowledge about the Milky Way, and in 
particular the discovery of substructure in the spheroid component, as we discussed above. 
Because of this, SDSS has been expanded to include a new project, the Sloan Extension for 
Galactic Understanding and Exploration (SEGUE). This extension will eventually include 
3500 sq. deg. of new imaging data which is collected in 2.5°-wide great circles on the sky, but 
these great circles do not in general match those laid out on the sky for the SDSS, and they 
are not adjacent to each other so they do not fill contiguous areas of the sky. The positions 
of these stripes were chosen to sparsely sample all directions in the sky that are visible 
from Apache Point Observatory, and include scans at constant Galactic longitude that pass 
through the Galactic plane. SEGUE will also obtain ~ 250, 000 spectra of Galactic stars to 
study stellar populations and kinematics. 

In this paper, we wish to use the photometry derived from imaging of Galactic stars to 
trace the three dimensional shape of the Galactic spheroid and the substructure, including 
tidal streams, contained within it. In this demonstration of the algorithm and the results 
from running it, we will focus on only one well-studied stripe of SDSS data: stripe 82 on 
the Celestial Equator. We use only F turnoff stars in 300 square degrees of this stripe, 
roughly centered on the place where the Sgr stream crosses it. Unlike galaxies, the distances 
of stars cannot be inferred from their radial velocities. By selecting only stars with colors 
consistent with spheroid F turnoff stars, we select a stellar population that is not a very 
good standard candle. But because we observe a large number of F turnoff stars we can 
find their underlying density distribution through statistical methods. We will assume the 
absolute magnitudes of the F turnoff stars in the population have a mean of M go = 4.2 and 
a dispersion of o~m jq = 0-6 . These are reasonable estimates for the Sgr dwarf tidal stream 



( iNewberg and Yannyil2006l ). 



1.3. The Technique of Maximum Likelihood 

Given a parameterized model and some data generated according to the model for some 
"true" values of the model parameters, the task of model estimation is to determine the set 
of parameters used in generating the data. Within a Bayesian setting, the model estimation 
problem can be reformulated as determining the a posteriori most likely parameters given 
the data and the model. In our case we have spatial positions for a set of stars in the 
spheroid, and a proposed model for the spheroid and a tidal stream which passes through 
it, and we would like to find the most likely values of the parameters in that model. 
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The likelihood of a set of parameters is the probability of obtaining a particular data 
set for a given set of parameters. Via Bayes' theorem, one can decompose the a posteriori 
probability of a particular set of parameters given the data and the model into the product 
of two terms, the first being the likelihood of the parameters and the second being the prior 
probability of the parameters. When the prior probability distribution over the parameters 
is uniform (as is typically assumed), the a posteriori probability is proportional to the 
likelihood, and so the maximum likel ihood technique may be used to find the most likely 



model parameters [see iFletcherl (119871 ) for more details]. 



There are two main tasks to implementing the maximum likelihood technique. The first 
is to develop the likelihood function, which measures the probability of observing a particular 
dataset given a parameterized model. The second is to maximize the likelihood with respect 
to the parameters. 

In developing the likelihood function, one typically assumes that each data point is 
independently generated, hence the likelihood of the dataset is the product of the likelihoods 
of each individual data point. Within our setting, we develop a parameterized model for 
the background Milky Way and Sagittarius stellar distributions from which we compute the 
likelihood. We prefer to maximize the logarithm of the likelihood, rather than the likelihood 
because the computation of the log-likelihood is numericall y more stable. We then use the 



standard technique of optimization via conjugate gradients (IFletcherl 119871 ) to maximize the 



likelihood and hence determine the most likely model parameters fitting the data. 

One advantage of the maximum likelihood framework is that the Hessian of the log- 
likelihood (at the maximum) gives the shape of the probability distribution over the pa- 
rameters, and hence allows us to determine the statistical error in the estimated model 
parameters. When the model parameters are physical quantities of interest, this statistical 
error translates to an error bar on the measured (i.e. estimated) physical quantity. 

In the remainder of this paper we present a new automated technique that can detect 
tidal debris in a spatial input catalog of stars from one stripe of SDSS data. Since the 
first stripe in which the Sgr dwarf tidal stream was detected in SDSS data was Stripe 82 



( INewberg. Yanny et al.ll2002l ). and it has therefore been the most well-studied; this is the 
section of data on which we test our maximum likelihood. Building upon this, we describe a 
technique to extract a catalog of stars that fits the density profile of the debris; this catalog 
can then be used to constrain the dynamical models. 
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2. Algorithm for Fitting Spheroid Substructure 

2.1. Goals of the Algorithm 

Starting with a parameterized density model for the spheroid, a parameterized density 
model for the Sgr dwarf tidal stream, the absolute magnitude distribution for F turnoff 
stars, and a dataset of observed (/, b) angular coordinates and apparent magnitudes for 
color-selected F turnoff stars in a 2.5° stripe through which a tidal stream passes, we find 
the parameters for a tidal stream and stellar spheroid concurrently that give us the highest 
likelihood of observing the data. This will allow us to tabulate summary statistics for the Sgr 
dwarf tidal stream, including position, width, and density. It will also allow us to test models 
for the spheroid component of the Milky Way, and ultimately, it will allow us to construct 
a catalogue of stars that fit the density profile of the Sgr tidal stream by probabilistically 
separating the stream from the spheroid. 

The parameters in our stellar density model describe the position, direction, and width 
of the tidal stream, as well as the shape of the smooth stellar spheroid. We start with an 
initial set of parameters and iterate using a conjugate gradient search technique until the 
optimum values for each parameter given the data set are determined. From the parameters 
which produce the maximum likelihood, we estimate the errors in each parameter. 

Once we have found the best parameters, we are able to use a separation algorithm to 
divide the data into two sets of stars, one of which has the density profile of the stream and 
the other of which has the density profile of the spheroid. While we cannot determine which 
of the individual stars in the data set belong to the stream and which belong to the spheroid, 
we can separate the input catalog into two catalogs that will have the density properties of 
the stream and the spheroid, respectively. These catalogs will allow for a close comparison 
of the stream density profile with that of simulations of tidal disruption and constrain the 
dynamical models used in simulations of tidal disruption. 

§2.2 gives a brief overview of the algorithm. §2.3 describes the construction of the prob- 
ability density function. §2.4 discusses the optimizaiton method we use. §2.5 describes the 
method by which we estimate the errors in the parameters. §2.6 discusses the application of 
a separation algorithm to distinguish stream and spheroid stars. §2.7 describes the compu- 
tational time required to run the algorithm and the application of parallel processing as a 
way to improve efficiency §2.8 discusses some limitations and possible enhancements that 
could be made in the future. A thorough discussion of the conjugate gradient and line search 
optimization method can be found in the Appendix. 
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2.2. Overview of the Algorithm 

The algorithm starts with initial guesses for the model parameters and an input catalog 
of observed F turnoff star positions as (I, b, g), which represent the Galactic longitude, in de- 
grees; Galactic latitude, in degrees; and reddening corrected q magnitude for each star. The 



appar ent magnitudes are corrected for reddening using the ISchlegel. Finkbeiner. fc Davis 



( 119981 ) reddening maps, as implemented in SDSS Data Release 6 (DR6). 



The essence of the maximum likelihood technique is to create a parameterized probabil- 
ity density function (PDF), and then find parameters in this function which maximize the 
likelihood of observing the data. We first select a parameterized stellar density function for 
the spheroid and its substructure. Since the parameterized models describe the density as a 
function of space (X, Y, Z), and the observations include only Galactic coordinates (I, b) and 
apparent magnitude, g, we need to relate spatial position to (I, b, g). The angular position in 
the sky is very well known for each star, but it is difficult to determine the distance to each 
star. We can only estimate the distance from its apparent magnitude, g. If all of the stars in 
the sample have the same absolute magnitude, then the conversion of apparent magnitude to 
distance is trivial. Unfortunately, this is not a good approximation for F turnoff stars, whose 
luminosity can deviate from the mean by a factor of 5 or more, even within the same age and 
metallicity population. Instead, we approximate the distribution of absolute magnitudes of 
F turnoff stars as Gaussian. 

The likelihood of observing the data is the product of the likelihoods of observing each 
of the stars within the dataset, given the model: 

N 

C(Q)=l[PDF(l u b l ,g l \Q), (1) 

i=i 

where the index i runs over the N stars, Q is a vector representing the parameters in the 
model, and the probability density function (PDF) is a normalized version of the stellar den- 
sity function that will be derived in the next section. Because the individual probabilities are 
small, we avoid numerical underflow by maximizing the average logarithm of the likelihood: 



1 1 N 

-ln£(Q) = -Y, ]nPDF &,bi,9i\Q), (2) 



which is maximized for the same parameters that the likelihood is maximized. 

We have chosen as input to this first version of the maximum likelihood algorithm a 
volume of the Galaxy sampled in a section of stripe 82. This volume is a piece of a wedge, 2.5° 
wide, with point at the Sun, limited in near and far end by the apparent magnitude range 
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16 < g < 22.5. For F turnoff stars with absolute magnitude M g = 4.2, the approximate 
distance range is 2.3 < R < 45.7 kpc from the Sun. The near and far edges of the volume 
are fuzzy due to the inexact relationship between apparent magnitude and distance, and 
the variable detection efficiency at the faint end. We chose this volume of a single stripe 
as the first application of the algorithm, for we do not have a global model for a general 
tidal stream; however, if we choose to look at a small enough volume of the stream (as is 
the intersection of the stream with Stripe 82), we are able to approximate its shape as a 
cylindrical density with cross-sectional Gaussian fall-off from its axis. We model the entire 
tidal stream as a set of these cylindrical stream segments. This model will not work if the 
stream is in the plane of the wedge, but we show later that this cylindrical approximation 
works for this stream in this wedge. 

The algorithm begins by calculating the average log-likelihood for an initial set of pa- 
rameters. A conjugate gradient with line search method is used to find the next set of 
parameters and the process iterates. Using this method, each subsequent set of parameters 
is guaranteed to produce an equal or higher likelihood in each iteration of the algorithm. The 
algorithm continues to generate new sets of parameters until convergence, that is defined by 
one of a number of conditions: either the gradient becomes so small that it would not pro- 
vide foreseeable improvement, the line search returns negligible movement, or a maximum 
number of iterations is reached. 



2.3. The Probability Density Function 

2.3.1. Tidal Stream Model 

In general, tidal streams follow a complex path through the sky. The stars in the 
structures may bunch up at apogalacticon and may have a complex cross-sectional density 
that varies with position along the stream. However, within a single SDSS stripe through 
which the stream passes it is reasonable to approximate the path of the stream as linear (see 
§3.3). 

We model the stream in a piecewise linear fashion with a separate set of parameters 
for each stripe of data in the SDSS, though in principle we could run this algorithm on any 
2.5°-wide great circle, or set of 2.5°-wide great circles on the sky. The length of the cylinder is 
limited by the edges of the data in one stripe. The cross section is circularly symmetric with 
a density that falls off as a Gaussian with distance from the stream axis. Figure 1 depicts 
the data volume and the relationship between the stream parameters and the segment of a 
cylinder that describes the stream. 
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The SDSS great circle coordinate system (//, v) is used to measure the angular position 
on the sky. v is bounded by —1.25° < v < 1.25° for all stripes, and measures the angular 
position across the narrow dimension of the stripe, v is by definition zero along the center 
of each stripe of data. \x measures the angular distance along the great circle swept out by 
each stripe. The inclination of each stripe is the maximum angle between that stripe and 
the Celestial Equator. Thus, fi, v, and the stripe inclination uniquely specify the angular sky 
coordinates. 

The vector c points from the Galactic center to the axis of the cylinder. Normally, a 
vector requires three parameters, one for each dimension. However, since we have the freedom 
to require that the position along the cylinder axis to which it points be in the plane that 
splits the volume in half along the narrow direction, we reduce that to two parameters. 
Enforcing the condition that v = 0, we are thus able to parameterize the stream center with 
only a radial distance from the Sun, R (in kpc), and the angular position along the stripe, 
/i (in degrees). Therefore, c(fj,,R) fixes the center point of a piece of tidal debris within an 
SDSS stripe and lies along the axis of our cylinder. 

The unit vector a describes the orientation of the axis of the cylinder. Again, we need 
only two parameters because this vector is constrained to be of unit length. We parameterize 
this vector using two angles 9 and 0, both in radians. 9 is the angle between a and the 
Galactic Z-axis. Z is perpendicular to the Galactic plane and points towards the North 
Galactic Cap. The azimuthal angle <fi is measured counter-clockwise around the Z 
viewed looking down on the Galaxy from the North Galactic Pole, starting from the X-axis, 
which points in the direction from the Sun to the Galactic Center. 

The last stream parameter, a (in kpc), specifies the stream width and is the standard 
deviation of the Gaussian distribution used to describe the density fall off with distance from 
the cylinder axis. For a star with spatial coordinates given by p, the distance, d, from the 
cylinder axis is given by 



In practice, this calculation is performed by first converting each vector to a Galactocentric 
Cartesian coordinate system. 

In summary, we use five parameters to define our cylindrical stream: //, R, 9, <p, and a. 
Figure 1 depicts how these parameters are defined with respect to the stripe volume. Using 
these parameters, the stellar density of the stream at point p is: 



Normalization of the stellar density will be considered once we have assembled the entire 
probability density function. 




(3) 



Pstream(p) OC e ^ . 



(4) 
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2.3.2. Spheroid 



We model the stellar spheroid with a standard Hernquist profile ( IHernquistl Il990l ) : 



Pspheroidip) oc — — - — , where (5) 
r{r + r ) 6 




and X, Y, and Z are Cartesian coordinates centered on the Galactic center. X and Y are 
in the Galactic plane, with X directed from the Sun to the Galactic center and Y in the 
direction of the Solar motion. Z is perpendicular to the Galactic plane and points in the 
direction of the North Galactic Cap. The spheroid is described by two parameters: ro which 
is a core radius, in kpc, and the dimensionless quantity, q, which is scaling factor in the Z 
coordinate direction. For q < 1 the spheroid is oblate, for q = 1 the spheroid is spherically 
symmetric, and for q > 1 the spheroid is prolate. 

There are many Galactic components (thick disk, bulge, etc) other than the spheroid 
that we could in principle add to our full model. However, we don't expect many disk stars 
in our sample because we are fitting data that is too far from the plane of the Milky Way, 
and the color-selected F turnoff stars are bluer than the turnoff of the thick disk stellar 
population. 



2.3.3. Absolute Magnitude Distribution 

In this section we address the fact that stars within our selected color range, which 
are primarily F turnoff stars, do not all have the same absolute magnitude. If we assume 
that all of the color-selected stars have the same absolute magnitude (equal to the mean 
absolute magnitude of the population) when we estimate their distances from the Sun, any 
substructure in the spheroid will appear to be elongated along our line of sight. To account 
for this, we calculate an "observed" spheroid spatial density that is elongated along our 
line of sight by convolution of the density model with the absolute magnitude distribution 
function along our line of sight. 

We model the distribution of absolute magnitudes as a Gaussian with a center of M go = 
4.2 and dispersion a go = 0.6, which is a si mplification of the F tu r noff s tar absolute magnitude 



distribution found for globular clusters by lNewberg and Yannyl (120061 ). Thus, letting M go be 
the absolute magnitude, 

M 9o = M go + AM go , (6) 
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where M go = 4.2 and AM So has a Gaussian distribution with zero mean and standard 
deviation 0.6. To account for this distribution over the absolute magnitude, we first derive 
the probability of observing a star per unit apparent magnitude per unit solid angle, assuming 
all stars are of absolute magnitude M go = 4.2, and then convolve in apparent magnitude with 
a Gaussian of dispersion 0.6 centered at zero. The result is the probability of observing a star 
per unit apparent magnitude per unit solid angle, with the absolute magnitude distribution 
taken into account. 

We will need to refer to several density functions in different spaces, and so we first 
lay down some definitions to make this discussion clear. We will refer to a density as pa{x) 
where A is the name of the density which will refer to the coordinate space in which it is 
defined and x refers to a generic variable in this space. Thus we define the 6 densities: 

px(x), p R (R,n), p 94 . 2 (94.2,ty, p go (g ,Q), p-ji(K(go),ty, p Xc (x), (7) 

where the first three densities are the model in Cartesian, spherical, and apparent magnitude 
coordinate systems; 

,_ x dV . . dV , . dV . . 

Px{x) = i^dz P *^w = im> ^ P ^ 9 "> n) = (8) 

The densities which take into account the distribution on the absolute magnitude space are 
denoted by p go (g ,fl), pfi{7l(go),Q), and p Xc (x), where the subscript c stands for convolved. 
All coordinate systems are centered at the Sun, the solar position is assumed to be 8.5 kpc 
from the Galactic center, and the direction from the Sun to the Galactic center is in the 
positive X direction. Here we have made a distinction between R and 1Z: R denotes the 
actual distance each star is from the Sun; TZ(go) denotes the distance a star of apparent 
magnitude g would be from the Sun if it had an absolute magnitude of M go = 4.2. We have 
also made the distinction between g and 54.2: go denotes the actual reddening corrected 
apparent magnitude of a star. ^4.2 denotes the apparent magnitude a star at distance R 
would be calculated to have if it had an absolute magnitude of M go = 4.2. We will adopt 
these definitions for the remainder of this paper. 

The Galactocentric Cartesian density, px, is the actual spatial density of stars as de- 
scribed in Equations 4 and 5 for the stream and spheroid, respectively. We need px c , the 
observed Galactocentric Cartesian density that is elongated along our line of sight to account 
for Gaussian distribution of absolute magnitudes with dispersion AM go . We will obtain p Xc 
through the sequence of transformations 

Px(x) -> p R {R,Q) -> p g42 {g 4 .2,ty -> p go (go,ty -> Pn(R>(go),n) -»• px c {x). (9) 

The relationship between these densities is determined by the transformations which take one 
coordinate space to the other. The X — > R mapping is the well known spherical coordinate 
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transform, 

p R {R,0) = R 2 p x {x). (10) 

If all of the stars had an absolute magnitude M go = 4.2, then one would measure an 
apparent magnitude 

g 4 . 2 = 4.2 + 51og 10 ( i g 5 ), therefore (11) 
R = n(g 4 . 2 ) = io°- 2 (^. 2 -4.2-io) (kpc), and 
dR = l -^Rdg 4 . 2 . 

Thus, the relationship between p R and p g is given by 

/ doAo , _ In 10 In 10^,o, s . . 

P 3 , 2 (^4. 2 ,fi) = -f^p R (R,n) = —R* px (x) = —<R?{g 4 . 2 ) Px {x). (12) 

The measured g is given by 

g = (74.2 + AM 90 . (13) 

Since g is the sum of independent random variables, its density is the convolution of the 
two densities, i.e., we have that p go (g ,Q) = p g * P AM go {go,ty, where the convolution is in 
the ^-dimension. Thus, 

/oo 
dg p 94 2 (g,^W(go- g;u), (14) 
-oo 

where Af is the Gaussian density function given by: 

N(x;u) = —^=e^, (15) 



with u = 0.6. Switching back from apparent magnitude to spherical coordinates, 

Since the coordinate spaces X c and 1Z are related by the spherical coordinate transformation, 
we may now collect our results and write the convolved density px c in terms of p x as follows: 

pxM = ^i^y M^(ffo),^), (17) 



K 3 (g )\nl0 

5 

^ 3 ((7o)lnlO 
1 Z" 00 



dg p g (g,ty ■ Af(g -g;u), 
c^ 3 (<?) • p x {x{n{g),Q)) ■ M(g - g;u), 
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where (H(go), fl) are the angular coordinates of x. 

We have derived the convolved stellar density function px c , for a generic stellar density 
Px-i which could represent either the stream or the spheroid densities in our present context. 
This convolution is performed separately on the stream and spheroid densities, to compute 
the functions p c s °ream(^b,7l(go)) a nd p^^J l, b, TZ(qn)). We us the numerical integration 
technique of Gaussian quadrature (IHeathl 120021 ) to do the convolution integral. 



2.3.4- Combined Probability Density Function 

We are now ready to compute the probability density function for the combination of 
the stream and spheroid densities. To do this we need the stellar densities of the stream 
and spheroid as derived in §2.3.3, the detection efficiency for finding stars as a function of 
apparent magnitude, the volume over which the density is defined, and one more parameter, 
e, which describes the fraction of stars in each of the two components (spheroid and stream). 

The detection efficiency function, £, wa s derived by fitting a sigmo id curve to the 



detection efficiency measurement in Figure 2 of lNewberg. Yanny et al.l (120021 ). The efficiency 



function accounts for the decrease in detection efficiency at faint magnitudes: 

£(qo) = — 7 — -i , where (18) 

s = (0.9402,1.6171,23.5877). 

The dimensionless parameter that defines the fraction of the input stars that are in the 
stream and the fraction that are in the spheroid is e. The parameter e is modeled, as with 
the stream parameters themselves, separately for each stripe of data that is analyzed. Thus, 
the value of e for a given stripe of data gives only the relative number of stars that comprise 
the stream as compared to the spheroid for that stripe of data and does not measure the 
fraction of stars in the stream as a function of position within the Galaxy. 

By definition, the fraction of stars in the stream is: 

e e 

P[stream] = -^—^ — -■ (19) 

Similarly for the spheroid 

P[spheroid] = :j— p — -■ (20) 

We introduce here the concept of constrained and unconstrained variables. Instead of e, we 
could have defined a variable / which is the fraction in the stream, and then the fraction 
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in the spheroid would be 1 — /. However, we would then need to maximize the likelihood 
subject to the constraint that this parameter must be between zero and one. To avoid that 
constraint, we instead introduced e, which is not constrained. If e is oo, then all of the stars 
are in the stream, if e = the stars are evenly split between the components, and if e = — oo 
then all of the stars are in the spheroid component. The other parameters are naturally 
unconstrained and do not need a conversion. 

In total, then, we fit eight model parameters: five from the stream fragment in the 
current stripe (p, R, 9, 0, a), two from the spheroid (r , q), and one that specifies the fraction 
of stars in each component (e). Together, we will refer to these parameters as the eight 
components of the parameter vector Q. 

The probability density function, then, is given by: 

£(K(g ))p™ am (l,b,K(g )\Q) 



PDF(l,b,K(g )\Q) 



1 + e*j£(K(g ))pr t ? eam (l,b,K(g )\Q)dV 



| 1 S(n(go))p? p n hermd (l,b,n(g )\Q) 

1 + # mn90))pTheroiA W<7o)|Q)^' 

where the integral is over the entire volume probed by the input data. 

Calculation of the volume integrals dominates the runtime of the algorithm. Since these 
integrals cannot be solved analytically for our functions, we define an integration mesh which 
divides the total volume into many volume elements. The edges of the volume elements are 
along constant g, p, or v. We then calculate the spheroid and stream probabilities at the 
center of each volume element and multiply this by the volume of that element. The values 
for the spheroid and stream probabilities in all volume elements are then summed and will 
be hereafter referred to as the stream and spheroid integrals, respectively. 



2.4. Optimization 

Now that we have defined a probability distribution function and a likelihood measure, 
we need to find the parameters that maximize that likelihood. At the beginning of the first 
iteration of the program, the likelihood is calculated for the input parameters. Thereafter, a 
new set of parameters and likelihood is calculated using a conjugate gradient search coupled 
with a line search technique. For a detailed description of the conjugate gradient and line 
search methods employed see the Appendix §A.l and §A.2, respectively. 
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2.5. Errors in the Parameter Estimates 

The accuracy of the parameter estimates depends upon the the shape of the likelihood 
surface at its maximum, which is governed by the number of stars in our input catalog. The 
fewer stars, the wider the peak at the maximum and the larger the statistical error. It is 
also limited by the accuracy with which we are able to numerically determine the maximum. 
We will address this latter point when we discuss the test data in §3. Here, we explain how 
we calculate the accuracy with which we expect to be able to determine each parameter. 

We assume that the likelihood surface for the parameters can be reasonably estimated 
by a Gaussian near the maximum. Then, an estimate for the variance of each parameter 
can be found from the square matrix of second order partial derivatives of the log-likelihood 
function, evaluated at the maximum. This matrix is called a Hessian Matrix, H. The 
variance matrix V is the inverse of the Hessian Matrix, normalized by the number of stars 
in the data set: 

V = —H 1 . (22) 

The square roots of the diagonal elements of the matrix provide us with the statistical error 
in the measurement of each parameter. We compute the Hessian numerically using a finite 
difference method (as we did the gradient) as follows: 

H}, — Hf, — H% + Hf, 
H i:i = ^, where, (23) 



H lj = C (Qj + h h Qi + ^)lall other Q k fixed > 

H ij = £(Qj ~ hj, Qi + /ii)l a ll other Q k fixed > 

H fj = £(Qj + h j, Qi - h i)\ a \\ other Q k fixed > 

H tj = C (Qj ~ h v Qi ~ h i) lall other Q k fixed ■ 

Here, Qj is the j th component of parameter vector Q, similarly for Qi, and hj is the pertur- 
bation value for Qj, similarly for hi. In practice we use the same h k to calculate the gradient 
as the Hessian. These h k can be seen in Table I. 



2.6. Separation 

Once we have found the values of the parameters in our model, we can use the model 
to separate the input catalog of stars into two catalogs, one having the density profile of the 
stream and the other having the density profile of the spheroid. It is possible to determine 
the probability that a star is in the stream or in the spheroid, however, it is not possible 
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to determine which stars are in the stream and which are in the spheroid. For example, 
suppose a star selected from the data set is computed to be a stream star with probability 
0.6 and that it is a spheroid star with probability 0.4. We would put that star in the stream 
catalog with probability 0.6 and in the spheroid catalog with probability 0.4, but there is a 
48% chance that the star was put in the wrong file. The separation has proved very useful 
for evaluating the effectiveness of the model fits to the data, and hopefully will be useful 
for fitting tidal stripping models to the stream density profiles. However, it may not be 
the best method for selecting spectroscopic follow-up targets to measure the composition 
and velocities of stream stars, for example. Spectroscopic follow-up targets might be better 
selected based on the probability that it is a stream star, rather than on the catalog to which 
it was randomly assigned. 

To perform this separation, we calculate for each star the probability, T, that it was 
drawn from the stream population given the parameters, where T is defined by 



spheroid portion of Equation 21. 

A uniform random number between and 1 is then generated and tested against this 
probability. If the probability of being drawn from the stream population is greater than 
that of the random number, then the star is placed into the stream star catalog; otherwise, 
it is placed within the spheroid catalog. We repeat this process, generating a new random 
number for each star until all stars have been placed into their appropriate catalog. Thus, 
we get two distinct populations: one that is a smooth spheroid and one that has all the 
density properties of the stream. 

We use this nondeterministic approach of testing the star probability against a random 
number for extracting the stream from the spheroid stars in order to get an accurate repre- 
sentation of the density profiles. Since the stream probability definition is based solely upon 
the star's distance from the stream axis, assigning all stars with a probability greater than 
a set value to the stream would result in assignment of all stars within a certain distance 
to the stream. This would not be an accurate representation of the stream as we would be 
carving out a cylinder of stars from the data rather than a set of stars that fit a specific 
distribution with Gaussian cross section in density. 





(24) 




where S(l,b,TZ(g ) | Q) is the stream portion of Equation 21, and B(l,b,lZ(g ) | Q) is the 
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2.7. Runtime and Distribution 

In determining the likelihood, a trade-off needs to be made in the accuracy of the 
numerical integral calculation in Equation 21. The smaller the step size, C, in each dimension 
the more accurate the integral calculation is, which can lead to faster convergence of the 
search method being used. However, increasing the precision of the integral calculation leads 
to polynomially longer calculation time (scaling like C 3 ). Thus, the precision to which the 
integral is calculated is the dominant factor in the runtime and accuracy of our astronomical 
model. 

Testing has shown that in order to generate a meaningful likelihood value in stripe 82, 
a minimum of 1.4 million points within the stream and spheroid integrals is required. On 
a single processor, this takes approximately 15 minutes, assuming there is no magnitude 
distribution and the non-convolved probabilities are used. As performing a single conjugate 
gradient descent can take thousands (or more) integral calculations, these long execution 
times presented an obstacle to this research. 

To generate a meaningful value for the magnitude convolution integral, each integral 
point requires the calculation of at least 30 surrounding points. This increases the calculation 
time of a single integral to 7.5 hours. Again, assuming a bare minimum of 1,000 likelihood 
evaluations for a single conjugate gradient descent (typically the range is closer to 5,000), 
evaluating one astronomical model would take almost a year on a single processor. 

To alleviate these massive computational requi rements, a Generic Maximum Likelihood 
Evaluator (GMLlH) has been developed and used (iDesell. et all 120071 ). GMLE allows for a 
likelihood calculation to use various distributed computing environments, such as multiple 
processors in a cluster, a heterogeneous grid of clusters, or even a supercomputer. 

Using this distributed framework with 88 processors on the Rensselaer Grid resulted 
in a 65 times speedup over a single processor and using 512 nodes of an IBM BlueGene/L 
system resulted in a 148 times speedup. It now takes only a few days to compute the 
maximum likelihood parameters for a particular model and stripe of data. In both the 
grid and supercomputing environments used, overhead of communication was very small 
compared to the calculation time. Communication time consists of approximately 1-10% 
of the time to calculate the likelihood and does not noticeably increase as more processors 
are used. With these observations, the calculation should scale to at least 1000 processors 
on a grid, and 10,000 processors on the IBM BlueGene/L. More precise calculations of the 
integral and magnitude convolution will allow the distributed calculation to scale to even 



1 GMLE is available for download along with more information at http://wcl.cs.rpi.edu/gmle 
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larger numbers of processors, due to the increased calculation time assuming relatively fixed 
communication times. 

GMLE has also been extended to allow maximum likelihood evaluation over the Berke- 
l ey Open Infrastructu re for Network Computing (BOINC) Internet computing framework 



( lAnderson et al.ll2005l ). This allows users to volunteer computing resources by download- 



ing a BOINC client and attaching to the MilkyWay@Home project^. Currently, over 500 
computers from around the world are being volunteered for the project and performance is 
comparable to 512 nodes of the IBM BlueGene/L. 



2.8. Limitations and Enhancements 

Generically, all enhancements to this algorithm either improve the Galaxy model, in- 
crease the quantity of data to which the model is fit, or improve the accuracy of speed and 
accuracy of the convergence. We plan to improve the algorithm in all of these aspects. 

We plan to build up a model of the Milky Way stellar density, starting with the spheroid 
component. Over time, we will be able to include additional streams and additional pieces 
of the Sagittarius dwarf tidal stream, a more sophisticated model of the smooth component 
of the spheroid, and disk components. The spheroid substructure is best fit in a piecewise 
manner as we have done in this paper, but the smooth components require that more direc- 
tions be fit simultaneously. To be able to fit multiple stripes during a single optimization 
it must be possible fit multiple debris pieces during a single run as well. To add additional 
tidal debris, we need only to add a term for each additional stream to the PDF: 



PDF 



h a 



strearrii i 



strearrii 

)dV 

1 Pslheroidi^ b, 7Z(go)\Q spheroid) 

1 + Etx *i) / PTheroiS, 6, K{g Q )\Q spher(M )dV' 
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where i and j denote the i th and j th stream of k total streams, respectively. We add five new 
parameters, Qstream.ii an d a sixth parameter, q, for each new stream segment. This change to 
the PDF would allow k pieces of tidal debris to be fit within a single or multiple stripes. We 
are considering whether we should enforce continuity conditions between adjacent sections 
of the tidal stream. Since it has been suggested that the smooth portion of the spheroid may 



2 http : / /milky way. cs . rpi . edu/ milky way / 
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not be symmetric, we would li ke also to increase the num ber of parameters in the smooth 



spheroid component, following iNewberg and Yannyl (120061 ) 



We are currently fitting only F turnoff stars, and assuming that the absolute magnitude 
distribution of these stars is a standard Gaussian for all populations. Ideally, we would use 
stars in a wider color range, and develop a probability distribution in magnitude and color for 
each population of stars. This would increase the number of stars available to the algorithm 
and would allow us to use additional information about real structures in the Milky Way - 
insisting that each population with a turnoff also has a main sequence and giant branch, for 
example. Implementing this requires significant development not only in the computation 
of the likelihood, but also in characterizing Galactic stellar distributions. Small steps in 
improving this aspect of the algorithm are currently underway. 

All of the preceding improvements to the algorithm increase the number of fit param- 
eters and the complexity of calculating the likelihood function. In parallel with scientific 
development of the algorithm, we will be improving the speed of convergence. We are ex- 
perimenting with other methods of finding the maximum, and increasing the number of 
architectures on which our software can be run in parallel. Future architectures that we 
will use include supercomputers and BOINC, which among other products supports the 
SETI@home application. 



3. Simulated Data 
3.1. Generation 

To test our algorithm we generated a simulated version of SDSS stripe 82. The "True 
Value" column of Table 2 shows the parameter values used for the data generation while 
Figure 2 shows a plot of the simulated data. The data set contained a total of 205,708 
simulated stars; 28,498 (13.85%)of which are stream stars. The stream model parameters 
and star ratio were chosen to be the same as the real data for stripe 82. The techniques of 
generation are described below. 



3.1.1. Stream 

We create the simulated stream using an active generation technique, which means 
that all points that are generated are valid stream stars. First, we define a set of stream 
parameters to generate over. From these parameters we then calculate the c and a vectors 
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defining our stream position and direction. We then generate three random numbers: a 
uniform random number that determines where along the stream axis the point is, and two 
Gaussian random numbers (with standard deviation a as defined in the parameters) to define 
the stream cross-sectional coordinates. The point is then converted to Galactic longitude, 
latitude and apparent magnitude. Thus far, we have assumed a constant absolute magnitude 
for each star in the stream. To account for the fact that there is an absolute magnitude 
distribution, we add or subtract a small amount from the apparent magnitudes of each of 
the stream stars in the dataset. A Gaussian random variable with standard deviation 0.6 
is generated (one for each point); this value is added to the apparent magnitude of the 
generated point. Finally, we keep points with a probability given by the efficiency function 
to produce a data set for a simulated stream. Stars near the magnitude limit of the data are 
tossed out of the sample with a probability equal to the efficiency of our object detection as 
a function of apparent magnitude in the real data. 



3.1.2. Spheroid 

The simulated spheroid is generated using a rejection sampling technique. Here we 
cannot apply the same active generation technique we used for the stream because the density 
distribution in the smooth component of the spheroid is more complex than a Gaussian. 
Instead, we define a rejection technique that allows us to generate over a given stripe. We 
cannot simply generate values for /i, u, and R, however, because the stars are not uniformly 
distributed in these variables. The distribution in v varies according to the cos(u), and the 
distribution in R must take into account the growing volume at increasing distance, ji may 
be generated uniformly over the stripe as it corresponds to a longitude, and the volume of 
space in each longitude bin is the same. The volume element in a stripe is given by: 

dV = R 2 cos(u)drdfidu, (26) 

and the corresponding stripe volume is: 

V = - »-)(sm(v + ) - sin(O), (27) 

where R max is the maximum radial distance, and the positive and negative superscripts refer 
to the maximum and minimum values of that coordinate for the stripe. These equations are 
used to define the functions which allow us to generate in the other two coordinates: 

v = sin _1 (u * (sin(^ + ) — sm(u~)) + sm(u~)), (28) 

R = RmaxU)3, 
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where u and iu are uniform random numbers between zero and one. \i is generated by a 
uniform random number between the stripe limits. Once a point is generated, its probability 
is calculated using the likelihood function developed in §2.3. The probability is then divided 
by the total probability possible for a star. The total probability possible for a star is simply 
the maximum value that can be returned by the likelihood function given the current volume 
and parameter set. The value of this ratio is then tested against a uniform random number 
between zero and one; the point is put into the data set if the value is greater than the value 
of the random number. 



3.1.3. Accuracy of Parameter Determination 

Using the method described in §2.5 we calculated the expected accuracy for each pa- 
rameter. Since we created the simulated data with a known set of parameters, we therefore 
know the true parameter values of the dataset. We can calculate the Hessian Matrix at the 
optimum values and subsequently get a set of error bars at the optimum values. These error 
bars therefore give us the accuracy with which we expect to find the values of the parameters. 
The values of these errors can be found in the column "Expected Deviation" (the statistical 
error bar) in Table 2 for the 205,708 star simulated data set. 

In addition to the statistical errors, there can be numerical errors driven primarily by 
the accuracy with which we determine the integrals. These numerical errors were estimated 
by determining the accuracy the likelihood can be calculated for a fixed number of points in 
the numerical integration and then calculating the 'error' in each parameter by perturbing 
a single parameter, holding all others fixed, until a change in the total likelihood is observed 
that is greater than the minimum accuracy threshold of the numerical integration. It should 
be noted that the values quoted as "Numerical Error" are heuristic errors, meaning they are 
estimated and not a true error bar. 



3.2. Testing 

We tested our algorithm by letting it optimize to convergence for eight randomized sets 
of input parameters. Randomized here means a random perturbation of the actual values 
by 75 percent of the parameter's value. Although parameters outside of this range are in 
principle allowed, the likelihood surface for parameter values very far from the correct values 
is so flat, and the gradient so small, that numerical errors dominate our measurement of 
the gradient. In these eight optimizations: five optimized to the true values, two optimized 
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to a local maxima which has a lower likelihood than that of the true values, and one did 
not converge because it was far enough from the correct parameter values that the gradient 
could not be accurately determined. The local maxima appears to be related to the fact 
that the stream crosses the stripe at low inclination, the algorithm sometimes fits the data 
with a much higher inclination angle and larger width to compensate. All other parameters 
are still optimized to their true values at this local maxima. 

The average values of the returned parameters from the five optimizations to the true 
values are presented in the "Optimized Value" column of Table 2, and the "Actual Devia- 
tion" is the deference between the "True Value" and the "Optimized Value." The returned 
parameters have very little deviation. The standard deviation of the five optimizations were 
calculated to be around an order of magnitude lower than the statistical error calculated 
from the Hessian. These can be be seen in the "Std. Deviations of Optimizations" column 
of Table 2. Note that if the theoretical deviation is large for a particular parameter then 
near the maximum, the likelihood changes little with the variation of that parameter. This 
leads to a larger "Actual deviation," and can also make it more difficult for the algorithm to 
numerically find the maximum since the gradient will be very small. The errors calculated 
from the Hessian assume that the maximum of the PDF is found exactly. 

As can be seen in Table 2, all of the parameters have "Actual Deviations" which are 
smaller than their "Theoretical Deviations" with the exception of tq. It has been discovered 
that the likelihood changes so little compared to perturbations in these parameters when 
close to the maximum that this parameter is not being calculated accurately enough to 
reach the "True Value" seen in Table 2. The true error bar should be taken as the sum, in 
quadrature, of the "Expected Deviation" and "Numerical Error" columns. Once this is done 
the "Optimized Value" of all parameters are comparable, within the errors, to the "True 
Value." 

Finally, we ran our separation algorithm upon the average returned parameter values 
to create two catalogs of stars given those parameters, and we plotted these separately in 
Figure 3. Clearly visible is the stream with a Gaussian density fall off in the left plot, while 
the right plot depicts a smooth Hernquist profile. 

We also tested our algorithm on this same simulated dataset of 205,708 stars without 
using an absolute magnitude correction, thereby simply assuming that all F turnoff stars 
have an absolute magnitude, M go = 4.2. This was done to see how important this correction 
was and to see whether the algorithm could be effectively used without it. The results showed 
that the parameters of /j, and r were the only values not affected by this correction; the rest 
of the parameters deviated wildly from their actual values with the stream parameters R, 9, 
and (f> being the worst of the eight parameters reaching upwards of 34a in error. 
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3.3. Robustness of Models 



The previous section showed that the algorithm produces the correct answer given that 
the data is drawn from our model. In this section, we test the algorithm on data that was not 
generated according to the assumed spatial model. We c onstr ucted a reproduction of the Sgr 
Dwarf Galaxy tidal disruption produced by lLaw et al.l (120051 ) using a semi-analytic N-body 
approach. Using their parameters for the Galactic potential and kinematic va lues, we con 



struc ted an orbit for Sagittarius using the NEMO Stellar Dynamics Toolbox (ITeuben. P.J. 



19951 ). Specifically, we used the case of a spherical dark matter halo, with a velocity disper- 
sion of 114 km/s. We used a Plummer Sphere with one million particles as our model for 
Sagittarius, The sphere was evolved along the orbit for 3.18 Gyr until it reached the present 
position of the Sgr dwarf spheroidal. The resulting plot of the disruption of Sa gittarius is 
shown in Figure 4. It is consistent with the results obtained by lLaw et al.l (120051 ). 



We then selected the volume corresponding to SDSS stripe 82 from the N-body simula- 
tion and fit it using our model Hernquist profile spheroid and maximum likelihood technique. 
Using the output model parameters found by our algorithm, we created a new simulated 
stream as described in §3.1.1. Finally, we compared the cross-sections of this simulated 
stream with that of a similar stream generated with our density model. 

A cross-section of these two streams 1 kpc thick and centered at the returned stream 
center are plotted, along with a histogram of the stars within the cross-section, in Figure 
5. Figure 5 (left panel) depicts the N-body simulation stream while Figure 5 (right panel) 
is that of our simulated stream. It can be seen that the N-body simulation is indeed not 
Gaussian. The drawback to testing with an N-body model is that there are no correct 
model parameters with which to compare; but the center, as shown in Figure 5 (left), and 
the direction, as shown in Figure 4, are reasonable. In order to assess the error in the 
determination of the stream center, one would have to define what is meant by the center of 
an asymmetric distribution, but the algorithm's choice seems reasonable. 

Next, we estimate the curvature of the Sgr tidal stream in stripe 82. Taking th e distance 
to the center of the stream in stripe 82 to be 29 kpc (INewberg. Yanny et al.ll2003l ) we calcu- 
late that the 2.5°-wide SDSS stripes would be 1.3 kpc thick at t his distance. The length of the 
strea m in stripe 82, taking an inclination of 30° from the stripe (IFreese. Gondolo. fc Newberg 
20051 ). is then 2.5 kpc. We estimate the radius of curvature of the Sgr tidal tail in the orbital 
pl ane at this point to be 18 kpc by fitting a circle to the two southern stripes with detections 



m 



Newberg. Yanny et al.l (120031 ) . The deviation from linear ca n then be found to be d = 0.06 



kpc a t the edge and is very small when compared to the 6 kpc (IFreese. Gondolo. fc Newberg 
20051 ) width of the stream or the 6.74 kpc width found in this paper. The linear approxima- 
tion is quite reasonable for the Sgr stream in stripe 82. 
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Finally, we tested our algorithm with differing spheroid models. Because it was easier 
to change the spheroid model we fit than to regenerate many sets of spheroid stars, we 
modified the exponential component of the model Hernquist profile the algorithm fits to the 
data. Three tests were performed with the exponent of the Hernquist profile set to 3.5, 2.5 
and 2.0 as compared to the value of 3.0 used to generate the data. The results from these 
tests can be found in Table 3. In all cases, the stream parameters and the epsilon parameter 
saw little to no change (all optimizations were within the errors calculated above). As 
expected, the spheroid parameter values were incorrect, but represent the best parameters 
for the model profile fit to the simulated data within our generated volume. 



4. Results From SDSS Stripe 82 

After validating the algorithm on the simulated data, we then applied the algorithm 
to the stripe 82 data as described in §2.2. Figure 6 shows a plot of this data where a 
piece of the Sgr tidal stream is clearly visible as a dense structure on left of the image at 
(a,So)« (31,21.5). 

We selected from the SDSS SkyServer Data Release 6 all of the sources in stripe 82 that 
had magnitudes 16 < g < 22.5, colors of 0.1 < (g — r) < 0.3 and (u — g) > 0.4, were 
identified as point sources, and were not EDGE or SATURATED. Since the selected stars 
are far from the Galactic plane, the reddening can reasonably be estimated from the total 
amount of dust in that direction in the sky. The subscript "0" indicates that we selected 
reddening corrected magnitudes from the SDSS databases. The cut in (u-g) is made to 
remove low-redshift QSOs. We limited the angular length along stripe 82 to 109°, with 
right ascension in the range 310° < a < 360°, < a < 59°. Since stripe 82 is centered on 
the Celestial Equator, the edges of the stripe are ±1.25° in declination and the angle along 
the stripe is given exactly by right ascension. A small section of data with right ascension 
323.2° < a < 323.6° and declination —1.0° < 5 < —0.7° was removed in order to avoid a 
globular cluster. The final dataset includes 115,907 stars. 

As was done for the test data, a number of randomized sets of initial parameters were 
used (having the same 75% perturbation) and the average of the parameters of those opti- 
mizations that converged to a minima was taken. In this case ten total optimizations were 
performed, five optimized to a minima, four to a local minima with lower likelihood than 
the true minima, and one was again a poor starting point allowing for no optimization. The 
Hessian method was then used to calculate the error bars for these parameters at the true 
minima. The errors, as well as the average parameter values, are tabulated in Table 4. There 
are additional, unknown systematic errors for real data that were not present for the model 
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data. The quoted errors assume that the model is an accurate representation of the stellar 
density function, including our formula for the smooth spheroid, tidal stream, and absolute 
magnitude distribution. 

The average parameter values were used as input for our separation algorithm to create 
stream and spheroid catalogs for stripe 82. These catalogs were then plotted in Figure 7. 
As can be seen, the Sgr tidal stream has been clearly extracted from the spheroid in the left 
panel while the spheroid remains in the right panel. 



We compared the average stream center and direction with Figure 4 from INewberg et al. 



( 120071 ) in Figure 8 of this paper. As can be seen in the figure, the stream direction is roughly 
tangent to the stream of 2MASS M giants (though at a different distance — there is a 
difference in scale between distances measured with the M giants and distances measured 
with F turnoff stars and RR Lyraes). Notice also that the stream direction does not lie exactly 
in the plane of the Sgr dwarf, but is tilted slightly with respect to the Xs gr — Ys g r plane (3° ± 



1°). W e also compare our direction to the directional results of iFreese. Gondolo. fc Newberg 
( 120051 ) . who estimated the angle between the observational plane and the Sagittarius dwarf 
stream to be 30° and the angle between the normal to the line of sight towards a point and 
the tangent of the stream at that point to be 10°. We calculated the former angle from 
our parameters by taking the angle between the observational plane of stripe 82 and the 
directional vector. Doing this we get a value of 30° ± 3°. Our value agrees perfectly. The 
second angle was calculated by first finding the line of sight vector to the center of the stream, 
calculating a normal to that, and then finding the angle between the stream direction and 
that normal. We found a value of 22.5° ±0.2°. This difference is explained by the drastically 
improved accuracy of the maximum likelihood algorithm over the estimate by eye for this 
more difficult angle. A tabulation of the stream center position and the stream direction in 
various coordinate systems can be found in Table 5. Note that the Cartesian Center and 
Cartesian Direction correspond to the vectors c and a, respectively, that are used in §2.3.1. 



Our results are in reasonable agreement with those of iNewberg. Yanny et al.l (120031 ). 
but with much tighter error bars. The center of the stream is shifted from a = 33.99° ± 1° to 
a = 31.3 7° ± 0.25°. This shift is slight ly larger than expected, but not overly so. Note also 
that the INewberg. Yanny et al.l (120021 ) estimate of the center in F turnoff stars is a = 33°, 
which is closer to our value. Our detection of r adial distance, however, corresponds exactly 
to that seen in the INewberg. Yanny et al.l (120031 ) . Our stream width of a = 2.86 is equivalent 
to a FWHM of 6.7 kpc, whi c h is i n good agreement with the value of 6 kpc estimated in 
Freese. Gondolo. fc Newbergj (120051 ). 



The stream ratio e provides a very good separation as shown in Figure 7 and cor- 
responds to a detection of approximately 16,000 stream stars out of the total 115,907 
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F turnoff stars in the strip e . By performing a calculation parallelling Equation 10 of 



Freese. Gondolo. &: Newbergi (120051 ) we are able to get an estimate of the total stellar mass 
in stripe 82 as a fraction of the Sgr dwa rf itself. We first calculate the number of stars in 
Figure 6 of iNewberg. Yannv et alT(l2002h with 0.1 < (g - r) < 0.3 and 16.0 < g < 22.5. 
This was calculated to be 2,194 stars . We then use this in place of the value calculated 
for F/G type stars in Equation 10 of [Freese. Gondolo. fc Newbergi (120051 ) to get the total 
number of F turnoff stars currently in the Sgr dwarf itself. This results in a total of 1,798,700 
F turnoff stars in the Sgr dwarf. Dividing this by the 16,050 stars found to be in the stream 
in stripe 82 yields that stripe 82 contains 0.9% of the Sgr dwarf's current F turnoff stars. 

Our spheroi d parameters of q = 0.46 and r = 19.4 are slightly unexpected given the 
previous work in iBell et al.l (120071 ) , which finds that the flattening parameter for the stellar 
halo closer to q = 0.6. However, this does not discount that the best likelihood for the stellar 
halo is obtained from these results, only that they were not the expected. Further results 
from other data sets and other stripes will be needed to determine the consistency of these 
results amongst other stripes and confirm such a strong flattening of the stellar halo. 



5. Conclusions 

We have measured the position, direction, distribution, width and number or F turnoff 
stars of the Sgr tidal debris in SDSS stripe 82 as well as the stellar spheroid parameters corre- 
sponding to flattening and core radius. Our work has resulted in the improved measurement 
of the Sgr stream width (FWHM = 6.74±0.06 kpc), position [{a, 5, R) = 31.37°±0.26, 0.0°, 
29.22 ± 0.20 kpc], and direction [{X, Y, Z) galactlc = -0.991 ± 0.007 kpc, -0.042 ± 0.033 kpc, 
0.127 ±0.046 kpc]. We measured the number of F turnoff stars in the stream in SDSS stripe 
82 (13.86 ± 0.06% of 115,907). There are 16,050 Sgr F turnoff stars in stripe 82 compared 
to 1.8 x 10 6 turnoff stars in the Sgr dwarf itself (0.9%). We also calculated values for the 
stellar spheroid parameters: the flattening parameter (q = 0.46 ± 0.024) and the core radius 
(r = 19.40 ± 0.59). While these spheroid parameters are not the expected values, they are 
intriguing and deserve further study. Finally, we were able to generate separate catalogs of 
stream stars and spheroid stars for the F turnoff stars in stripe 82. These catalogs, while 
not a true representation of the actual stars in the stream or spheroid, provide the correct 
density characteristics of these structures and can be compared with the density profiles of 
N-body simulations. 

The results were calculated using a new method to study tidal debris that is more robust 
than previous methods. This method will prove valuable to the study of tidal disruption and 
Galactic structure as it provides a means to accurately and efficiently study many variables 
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used to model the spheroid and tidal debris. Also, we are able to extract a catalog of stream 
stars which can be used to better constrain models of the Sgr dwarf disruption by providing 
actual stream stellar positions to compare with the data. 

The algorithm itself is quite flexible. The stream and spheroid models can be easily 
changed should more accurate ones be discovered. The new model definitions would simply 
replace the existing probability code with no other changes. Also, with minor changes to the 
volume definitions it would be possible to use data from surveys other than the SDSS as an 
input catalog. With some revisions, the code could be adapted to search for multiple debris 
and use multiple volumes to provide more accurate results. 
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APPENDIX 

A. Optimization Technique 

A.l. Conjugate Gradient 

The conjugate gradient method is similar to the gradient method in that we perturb 
the current parameter value, Qi, by a small amount, hi, and get a direction for the change in 
the parameter space. However, unlike the gradient method, the conjugate gradient method 
enforces the condition that each direction is conjugate to the previous direction. This means 
that as we move along the new search direction, the component of the gradient parallel to 
the previous search direction must remain zero. This condition speeds up the optimization 
by using search directions that are non-interfering. A further discussion may be found in 



Fletcher! (Il987f ). 



The i th component of the gradient vector, G, is calculated numerically by perturbing 
Qi, holding all other Qj fixed: 



_ C(Qi + hi) - C(Qi - hi 

Lri — 



2h, 



(Al) 

all other Qj fixed 



where Gi is the i th component of the gradient vector, Qi is the i th parameter, hi is the 
perturbation amount for that parameter, and C is the likelihood function from §2.2. 

This gradient is used to calculate the direction, D, for the first iteration: 

Di = Gi- (A2) 

We take the positive value of the gradient, for we want to maximize the likelihood. Should we 
want to perform a simple gradient ascent we would do this same calculation for every iteration 
and use solely the direction calculated above; however, we wish to apply the conjugate 
gradient technique which always calculates a direction conjugate to the previous directions. 
This is done by calculating a multiplier, B i} based upon the current gradient and the previous 
gradient, as; 

Bi = (A3 ) 

where Gi denotes the gradient vector for the i th iteration. This value is then used to calculate 
the new conjugate gradient direction: 

Di = Gi + Bi* Di-i, (A4) 
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where A denotes the direction vector of the i th iteration (the direction to move the current 
parameters) and Gi is the gradient of the current parameters, Q iy at the the i th iteration. 

Due to the large variety of parameters that we utilize, all on different scales, we do not 
use one value for our perturbation value, h, for it could potentially produce poor estimates 
for the gradient if the perturbation is too large or small for a given parameter. If h is too 
large, then the computed gradient is not an accurate approximation of the true gradient 
at the current point; if h is too small then we may encounter computational errors. We 
therefore choose an appropriate value for each parameter independently. Table 1 shows the 
perturbation values for each parameter. These numbers were chosen in two ways. The 
values of the gradient were calculated at a number of differing perturbation values in order 
to determine how much the gradient changed. Also, plots of all parameter spaces were 
generated in order to determine the distribution of each parameter. Combining these two 
techniques we were able to determine a value for the perturbation for each parameter that 
would be small enough to provide a very accurate gradient calculation, but also large enough 
to overcome any anomalous behavior within the parameter's distribution. 



A. 2. Line Search 

The line search technique seeks to determine the value, a, that minimizes the function 

$(a) = £(Q k + aD k ), (A5) 

where D k is the direction and Q k is the current set of parameters. After finding a* which 
minimizes we then update 

4+i < — Qk + a*D k . (A6) 

A bracketing method must be employed first to ensure that the minimum of the function 
is within the range we are searching. This is done by calculating the likelihood at three 
points along the current search direction: the current parameter values plus zero, one, and 
two times D k . This corresponds to the current position, the current position plus the full 
direction, and the current position plus twice the full direction. If the middle point does not 
have a likelihood greater than the two endpoints, the endpoint with the highest probability 
becomes the new center point, the center point the end, and a new third point is calculated by 
expanding to two times the current factor times direction and the new likelihood is calculated. 
Iteration continues until the midpoint has a greater likelihood than both endpoints. 

These three points are then passed to the line search algorithm which will use them 
to find the peak of the parabolic function defined by moving along the current directional 
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vector and calculating the likelihood at points along this path. The line search iterates as 
follows: first a guess is made for the value of a by fitting a parabola via the calculation: 

u = aj*(L 2 - L 3 ) + al*(L 3 - L^ + aj*^- L 2 ) : (A7) 

b = ai * (L 2 - L 3 ) + a 2 * (L 3 - Li) + a 3 * (Li - L 2 ), 
u 

a = 0.5*—, 
b 

where a±, a 2 , and a 3 are the factors multiplied by the current direction that produce the 
parameters with likelihood L ± , L 2 , and L 3 , respectively. For the first iteration, a±, a 2 , and 
a 3 are the values returned from the bracketing method, a is then the current guess for the 
value to minimize Equation A5. The likelihood L a is then calculated using the parameters 
generated using a. The new order of the three points is: 

if : a > a 2 and L a > L 2 ; then : 1, 2, a, (A8) 
if : a > a 2 and L a < L 2 ; then : 2, a, 1, 
if : a < a 2 and L a > L 2 ; then : a, 2, 3, 
if : a < a 2 and L a < L 2 ; then : 1, a, 2. 

In short, the distance along the search direction is reduced based upon the value of a and 
its corresponding likelihood. These new points are then used to calculate a new guess for 
a using Equations A6 through A8. Iteration continues for a set number of iterations until 
returning the optimum value of a; we currently use three. 

Once we have calculated the change in the parameters, utilizing the conjugate gradient 
and line search methods, we then update the value of the probability for the fit using the 
same method as discussed in §2.3 and continue to maximize our probability through the use 
of the conjugate gradient and line search methods. The algorithm stops iterating and returns 
the current values when the value of the largest component of the conjugate gradient drop 
below a set threshold, the line search returns a value that would cause negligible movement, 
or the algorithm has completed a maximum number of iterations. 
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Fig. 1. — Stripe and stream parameter definitions. The segment of a stream which passes 
through an SDSS stripe is cylindrical within an individual stripe, with density that falls off 
as a Gaussian with distance from its axis. The coordinates //, is, and r are used to define 
SDSS stripes. We adopt these coordinates to define a vector, c(fi,R) [y = 0), which points 
to the center of the stream from the Galactic center. We then define a stream directional 
vector a(9, (f>) of unit length. Finally, we define the stream width as a, which is the standard 
deviation of the Gaussian that defines the density fall-off of the stream. 

Fig. 2. — Simulated data wedge plot. We show a Sun-centered density plot of the 205,708 
stars generated to mimic SDSS stripe 82. g is labeled along the radial spokes; the circle 
indicates where go = 23. The angle fi in degrees is marked along this circle. For stripe 82 
H = a and v — 5. The simulated stream is easily discerible at (a, go) = (33.4, 21.4). 

Fig. 3. — Simulated stream and spheroid stars after separation. Here we have plotted the 
same 205,708 star simulated data set similarly to Figure 2. However, we have separately 
plotted the stream (left) and spheroid (right) catalogs returned from running our separation 
algorithm. The spheroid is recovered as a smooth Hernquist function after the removal of 
the stream component. 

Fig. 4. — Plot of the disruption of the one million particle simulated Sagittarius Dwarf with 
a spherical dark matter halo (q = 1), and a velocity dispersion of Vhaio = 114 km/s in the 
Sgr XY orbital plane. To reduce crowding, the stream was sampled one in one hundred. The 
thin line is the future orbit of the core, while the dotted line is the past orbit. The past orbit 
is not closed, but in fact follows the lower trailing stream. The arrow depicts the center and 
direction returned by our algorithm optimization. 

Fig. 5. — Cross-sections of Sgr N-body simulation stream (left panel) and simulated stream 
(right panel). The simulated stream matches our stream model and was generated from the 
parameters fit to the N-body simulation. The upper figure in both panels show the 1 kpc 
thick cross-section of the respective data set. The cross-sections are centered at the best-fit 
value of the optimization of the center in the N-body simulation from Figure 4. Here the axes 
are perpendicular to the stream direction. The X-axis is 0.053X + 0.055Y + 0.997Z and the 
Y-axis is 0.055X + 0.997Y — 0.058Z, where X,Y,Z are Galactocentric Cartesian coordinates 
with the Sun at X = -8.5 kpc and moving in the direction of positive Y. The lower figure of 
both panels is a histogram of those stars within the cross-section binned along the X-axis. 
The heavy dashed line shows a Gaussian distribution with standard deviation given by a 
from the fit to the N-body simulation. Note that the cross-section of the N-body simulation 
is non-Gaussian and somewhat asymmetric. Also note the model simulated stream is well 
fit by a Gaussian. The density distributions in the left and right panels are not the same, 
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however, the algorithm still fits a reasonable center for the non-Gaussian N-body data set. 

Fig. 6. — F turnoff stars within SDSS stripe 82. Here we plot, as in the Figure 2, F turnoff 
stars with color cut 0.1 < (g — r)o > 0.3 and (u — g)$ > 0.4 that are within the volume 
defined by 310 < a < 419, ±1.25 = 5, and 16 < g < 22.5 that were detected in SDSS stripe 
82. Sgr tidal debris is visible around (a, go) = (31, 21.5). 

Fig. 7. — F turnoff stars within SDSS stripe 82 after separation. Here we have plotted the 
data from Figure 6 after being separated into two distinct catalogs: one for stream and one 
for spheroid. The stream (left) has clearly been extracted from the spheroid (right). 

Fig. 8. — Plot of the stream center and dir ection in the Sgr dwarf plane compared to other Sgr 



detections. Here we reproduce Figure 4 of iNewberg et al.l (120071 ) which displays detections of 



the Sgr stream in A-colored stars in eleven stripes, where the filled circles and larger squares 
represent leading debris, the open circles trailing debris, and the smaller squares Sgr debris 
on the opposing side of the Sgr orbi tal plane. These are p lotted along with the positions of 



2MASS M giants from Figure 11 of iMajewski et al.l (120031 ) . The arrow shows our improved 



measurement of the position and direction in one place on the Sgr tidal stream. The length 
of the directional vector is arbitrary representing only the projection of an elongated unit 
vector onto the respective planes. Note that the direction is tangent to the 2MASS M star 
data, and plausibly on a smooth path from stripe 86 (the open circle just to the right of our 
detection for stripe 82) to stripe 27 (the open circle farthest to the left in the lower figure). 
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Table 1. Perturbation values used for the gradient and Hessian 



Parameter 




h 




H (dcg) 


3 


10" 


5 


R (kpc) 


4 


10" 


5 


6 (rad) 


6 


10" 


5 


<f> (rad) 


4 


10" 


5 


a (kpc) 


4 


10" 


6 


£ 


1 


10" 


6 




4 


10" 


6 


ro (kpc) 


8 


10" 


4 
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Table 2. Results for the simulated data set of 205,708 stars 



Parameter 


True 


Expected 


Numerical 


Optimized 


Actual 


Std. Deviation 




Value 


Deviation 


Error 


Value 


Deviation 


of Optimizations 


H (deg) 


31.361 


0.233 


0.05 


31.443 


0.082 


0.064 


R (kpc) 


29.228 


0.167 


0.04 


29.217 


0.011 


0.010 


6 (rad) 


1.445 


0.032 


0.003 


1.421 


0.024 


0.0005 


<j> (rad) 


3.186 


0.049 


0.001 


3.182 


0.004 


0.002 


(j (kpc) 


2.854 


0.033 


0.015 


2.858 


0.004 


0.009 


e 


-1.828 1 


0.005 


0.002 


-1.833 


0.005 


0.005 


<? 


0.670 


0.013 


0.000 


0.671 


0.001 


0.0004 


r (kpc) 


13.500 


0.276 


0.15 


13.917 


0.417 


0.016 



Corresponds to a stream star ratio of 13.85 percent: 28,498 in stream; 177,210 in spheroid. 
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Table 3. Stream fitting results for tests with incorrect spheroid model. 



Parameter 


True 


Expected 


a = 3.5 


a = 2.5 


a = 2.0 




Value 


Deviation 


Optimization 


Optimization 


Optimization 


H (deg) 


31.361 


0.233 


31.449 


31.526 


31.457 


R (kpc) 


29.228 


0.167 


29.094 


29.326 


29.108 


6 (rad) 


1.445 


0.032 


1.426 


1.452 


1.437 


<f> (rad) 


3.186 


0.049 


3.172 


3.168 


3.160 


a (kpc) 


2.854 


0.033 


2.869 


2.848 


2.865 
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Table 4. Results for the 115,907 F turnoff stars in SDSS Stripe 82 



Parameter 


Optimized 


Statistical 


Numerical 


Std. Deviation 




Value 


Error 


Error 


of Optimizations 


fi (deg) 


31.373 


0.244 


0.08 


0.008 


R (kpc) 


29.218 


0.184 


0.07 


0.012 


9 (rad) 


1.444 


0.044 


0.01 


0.001 


<fi (rad) 


3.184 


0.034 


0.008 


0.002 


cr (kpc) 


2.862 


0.025 


0.008 


0.009 


e 


-1.827 1 


0.005 


0.001 


0.000 


<? 


0.458 


0.023 


0.005 


0.001 


r (kpc) 


19.404 


0.581 


0.09 


0.051 



Corresponds to a stream star ratio of 13.86 percent of the 115,907 F 
turnoff stars; about 16,000 stream stars. 



-42- 



Table 5. Table of the stream center and direction as calculated for different coordinate 
systems. All distances are in kpc and angles are in degrees. 



Physical 


Coordinate 


Value 1 


Value 2 


Value 3 


Quantity 


system 








Center 


Equatorial (a, 8, R) 


31.37 


0.0 


29.22 1 


Center 


Equatorial (a, <S, go) 


31.36 


0.0 


21.53 


Center 


Galactic (I, b, R) 


159.223 


-57.558 


29.22 1 


Center 


Galactic (I, b, go) 


159.223 


-57.558 


21.53 


Center 


Cartesian (X, Y, Z) 


-23.154 


5.560 


-24.658 


Direction 


Cartesian (X, Y, Z) 


-0.991 


-0.042 


0.127 


1 Assuming 


an absolute magnitude of M 30 = 


4.2. 
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